!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! This module provides types and procedures to work with
!! permutations.
!!
!! \n
!!
!! A permutation is represented by the following fields:
!! \li\c d    : dimension
!! \li\c p    : parity \f$\pm 1\f$ (or \f$0\f$ for uninitialized)
!! \li\c i(d) : a permutation of the integer vector \f$[1,\ldots,d]\f$
!!
!! Permutations of dimension \f$d\f$ are ordered according to the
!! lexicographic order; this ordering is absolute. Given a permutation
!! \c pi, its index can be obtained as \c idx(pi).
!!
!! In addition to permutations, we also define "data permutations",
!! which essentially attach some arbitrary integer values to a
!! permutation. Data permutations are described as:
!! \li\c pi             : a permutation of dimension \c d
!! \li<tt>x(pi\%d)</tt> : integer values (the data)
!!
!! The field \c pi\%i of a data permutation can be left
!! undefined, in which case we have \c pi\%p.eq.0. Data permutations
!! can be ordered according to the lexicographic order, but this
!! ordering is not absolute, since there is no bound on the values
!! that can appear among the data \c x. For this reason, data
!! permutations don't have an index, but can only be compared as \c
!! dpi1.gt.dpi2, and similar operations for <tt>.ge.</tt>,
!! <tt>.lt.</tt>, <tt>.le.</tt> and <tt>.eq.</tt>. Notice that these
!! comparisons are based on the values \c dpi1\%x and \c dpi2\%x, and
!! not on \c dpi1\%pi\%i and \c dpi2\%pi\%i.
!<----------------------------------------------------------------------
module mod_perms

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_linal, only: &
   mod_linal_initialized, &
   sort, fsort

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_perms_constructor, &
   mod_perms_destructor,  &
   mod_perms_initialized, &
   t_perm, t_dperm, &
   operator(.eq.),  &
   operator(.ne.),  &
   operator(.lt.),  &
   operator(.le.),  &
   operator(.gt.),  &
   operator(.ge.),  &
   operator(*),     &
   operator(**),    &
   perm_table,      &
   idx,             &
   dperm_reduce,    &
   fact,            &
   comb_table

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !< permutation type
 !< \bug the field d should be turned into a LEN parameter
 type t_perm
   integer :: d     !< number of elements
   integer :: p = 0 !< parity \f$\pm 1\f$, or 0 for undefined
   integer, allocatable :: i(:) !< permutation of \f$[1,\ldots,d]\f$
 end type t_perm

 !< data permutation type
 !< \bug the field d should be turned into a LEN parameter
 type t_dperm
   type(t_perm) :: pi !< permutation
   integer, allocatable :: x(:) !< data (pi\%d)
 end type t_dperm

! Module variables

 ! public members
 logical, protected ::               &
   mod_perms_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_perms'

 interface operator(.eq.)
   module procedure eq
 end interface
 interface operator(.ne.)
   module procedure ne
 end interface
 interface operator(.lt.)
   module procedure lt
 end interface
 interface operator(.le.)
   module procedure le
 end interface
 interface operator(.gt.)
   module procedure gt
 end interface
 interface operator(.ge.)
   module procedure ge
 end interface
 interface operator(*)
   module procedure comp
 end interface
 interface operator(**)
   module procedure pow
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_perms_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_linal_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_perms_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_perms_initialized = .true.
 end subroutine mod_perms_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_perms_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_perms_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_perms_initialized = .false.
 end subroutine mod_perms_destructor

!-----------------------------------------------------------------------
 
 !> Compute all the \f$d!\f$ permutations of dimension \f$d\f$ in
 !! lexicographic order.
 !<
 pure subroutine perm_table(ptab,d)
  integer, intent(in) :: d !< dimension
  type(t_perm), allocatable, intent(out) :: ptab(:) !< permutation table

  integer :: df, i
  integer, allocatable :: p(:), pi(:,:)

   df = fact(d)

   ! build temporary arrays
   allocate(p(df),pi(d,df))
   call genperms(p,pi,(/ (i, i=1,d) /))

   ! copy the permutation in ptab
   allocate(ptab(df))
   do i=1,df
     ptab(i)%d = d
     ptab(i)%p = p(i)
     allocate(ptab(i)%i(d)); ptab(i)%i = pi(:,i)
   enddo

   deallocate(p,pi)
 end subroutine perm_table
 
!-----------------------------------------------------------------------
 
 !> Generate array permutations
 pure recursive subroutine genperms(p,pi,v)
  integer, intent(in) :: v(:) !< vector to permute
  integer, intent(out) :: p(:), pi(:,:) !< parity and perms.
  
  integer :: d, fdm1, i, cjs, cje
  integer :: skip_i(size(v)-1)

   d = size(v)

   if(d.eq.1) then
     p = (/ 1 /)
     pi(1,1) = v(1)
   elseif(d.eq.2) then
     p = (/ 1 , -1 /)
     pi(1,1) = v(1); pi(2,1) = v(2)
     pi(1,2) = v(2); pi(2,2) = v(1)
   else
     fdm1 = fact(d-1)
     skip_i = (/ (i, i=2,d) /)
     do i=1,d
       cjs = (i-1)*fdm1 + 1 ! start column index
       cje =  i   *fdm1     ! end column index
       pi(1,cjs:cje) = v(i)
       call genperms( p(cjs:cje) , pi(2:d,cjs:cje) , v(skip_i) )
       ! correct parity for even i
       if(even(i)) p(cjs:cje) = -p(cjs:cje)
       ! update index vector
       if(i.lt.d) skip_i(i) = i
     enddo
   endif
 end subroutine genperms
 
!-----------------------------------------------------------------------
 
 !> Get the absolute index of a permutation (using the lexicographic
 !! order)
 !<
 elemental function idx(pi) result(i)
  type(t_perm), intent(in) :: pi !< permutation
  integer :: i !< permutation index

  integer :: j,k,n

   ! Factorial counting: for "digit" j, we build smaller permutation by
   ! 1) fix j-th digit so that
   !  a) it's smaller than i(j)
   !  b) it hasn't been used yet (this is the inner loop on k)
   ! 2) fix the remaining d-j digits in arbitrary order
   i = 1
   do j=1,pi%d
     n = pi%i(j)-1
     do k=1,j-1
       if(pi%i(k).lt.pi%i(j)) n = n-1
     enddo
     i = i + n*fact(pi%d-j)
   enddo
 
 end function idx
 
!-----------------------------------------------------------------------

 !> Given an integer array \c v, sort it an return the result in a
 !! t_dperm object. The \c pi field contains the permutation used to
 !! sort the vector.
 !<
 pure function dperm_reduce(v) result(dperm)
  integer, intent(in) :: v(:)
  type(t_dperm) :: dperm

  integer :: i

  dperm%pi%d = size(v)
  allocate(dperm%pi%i(dperm%pi%d))
  allocate(dperm%x   (dperm%pi%d))

  dperm%pi%i = (/ (i, i=1,size(v)) /)
  dperm%x    = v
  call sort(dperm%x,dperm%pi%i,dperm%pi%p)

 end function dperm_reduce
 
!-----------------------------------------------------------------------
 
 elemental function ne(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: ne

   ne = any(dpi1%x.ne.dpi2%x)
 
 end function ne
 
!-----------------------------------------------------------------------
 
 elemental function eq(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: eq

   eq = .not.(dpi1.ne.dpi2)
 
 end function eq
 
!-----------------------------------------------------------------------

 elemental function lt(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: lt

  integer :: i

   lt = .false. ! to account for the case of equality
   check: do i=1,dpi1%pi%d
     if(dpi1%x(i).lt.dpi2%x(i)) then
       lt = .true.
       exit check
     elseif(dpi1%x(i).gt.dpi2%x(i)) then
       ! lt already set to .false.
       exit check
     endif
   enddo check

 end function lt
 
!-----------------------------------------------------------------------

 elemental function le(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: le

   le = (dpi1.lt.dpi2).or.(dpi1.eq.dpi2)
 
 end function le
 
!-----------------------------------------------------------------------

 elemental function gt(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: gt

   gt = .not.(dpi1.le.dpi2)
 
 end function gt

!-----------------------------------------------------------------------

 elemental function ge(dpi1,dpi2)
  type(t_dperm), intent(in) :: dpi1, dpi2
  logical :: ge

   ge = .not.(dpi1.lt.dpi2)
 
 end function ge

!-----------------------------------------------------------------------
 
 elemental function comp(pi1,pi2)
  type(t_perm), intent(in) :: pi1, pi2
  type(t_perm) :: comp

   comp%d = pi1%d
   comp%p = pi1%p*pi2%p
   allocate(comp%i(comp%d)); comp%i = pi2%i(pi1%i)
 
 end function comp
 
!-----------------------------------------------------------------------

 elemental function inv(pi)
  type(t_perm), intent(in) :: pi
  type(t_perm) :: inv

  integer :: i

   inv%d = pi%d
   inv%p = pi%p
   allocate(inv%i(inv%d)); inv%i(pi%i) = (/ (i, i=1,inv%d) /)
 
 end function inv
 
!-----------------------------------------------------------------------

 elemental function pow(pi,n)
  type(t_perm), intent(in) :: pi
  integer, intent(in) :: n
  type(t_perm) :: pow

  integer :: i
  type(t_perm) :: base

   if(n.lt.0) then
     base = inv(pi)
   else
     base = pi
   endif

   pow = id(pi%d)
   do i=1,abs(n)
     pow = base*pow
   enddo
     
 end function pow
 
!-----------------------------------------------------------------------

 pure function id(d)
  integer, intent(in) :: d
  type(t_perm) :: id

  integer :: i

   id%d = d
   id%p = 1
   allocate(id%i(id%d)); id%i = (/ (i, i=1,id%d) /)
 
 end function id
 
!-----------------------------------------------------------------------
 
 elemental function even(n) result(l)
  integer, intent(in) :: n
  logical :: l

   l = (2*(n/2)) .eq. n
 
 end function even
 
!-----------------------------------------------------------------------
 
 !> Compute \f$n!\f$
 elemental function fact(n) result(fn)
  integer, intent(in) :: n
  integer :: fn

  integer :: i
 
   ! max(1,n) is introduced to deal with n=0
   fn = product( (/ (i, i=1,max(1,n)) /) );
 end function fact

!-----------------------------------------------------------------------
 
 !> Compute all the \f$\frac{d!}{k!(d-k)!}\f$ combinations with
 !! \f$k\f$ elements of the \f$d\f$ elements of \c elements. The
 !! combinations are returned in lexicographic order.
 !<
 pure subroutine comb_table(ctab,elements,k)
  integer, intent(in) :: elements(:) !< elements
  integer, intent(in) :: k !< number of elements for each combination
  integer, allocatable, intent(out) :: ctab(:,:)!< combinations

  integer :: d, nc, sorted(size(elements))

   d = size(elements)
   if((k.lt.1).or.(d.lt.k)) then
     allocate(ctab(0,0))
   else
     nc = fact(d)/(fact(k)*fact(d-k))
     allocate(ctab(k,nc))
     sorted = fsort(elements)
     call gencombs(ctab,sorted,k)
   endif

 end subroutine comb_table
 
!-----------------------------------------------------------------------
 
 !> Generate array combinations
 pure recursive subroutine gencombs(combs,elements,k)
  integer, intent(in) :: elements(:)
  integer, intent(in) :: k
  integer, intent(out) :: combs(:,:)

  integer :: d, i, cjs, cje

   d = size(elements)
   if((k.lt.1).or.(d.lt.k)) return

   if(k.eq.1) then
     combs(1,:) = elements
   else
     cjs = 1
     do i=1,d-k+1
       cje = (cjs-1) + fact(d-i)/(fact(k-1)*fact((d-i)-(k-1)))
       combs(1,cjs:cje) = elements(i)
       call gencombs(combs(2:k,cjs:cje),elements(i+1:d),k-1)
       cjs = cje + 1
     enddo
   endif
   
 end subroutine gencombs

!-----------------------------------------------------------------------

 !> This subroutine tests the main functionalities of this module. It
 !! should not be used in real application, but it can be useful to
 !! check changes to the module. To use it, add it to the \c PUBLIC
 !! list and simply call it in a program like the following \code
 !! program test
 !!  use mod_messages
 !!  use mod_kinds
 !!  use mod_linal
 !!  use mod_perms
 !!   call mod_messages_constructor()
 !!   call mod_kinds_constructor()
 !!   call mod_linal_constructor()
 !!   call mod_perms_constructor()
 !!   call test_module()
 !!   call mod_perms_destructor()
 !!   call mod_linal_destructor()
 !!   call mod_kinds_destructor()
 !!   call mod_messages_destructor()
 !! end program test
 !! \endcode
 !<
 subroutine test_module()

  integer :: i, d
  integer, allocatable :: v(:), combs(:,:)
  real, allocatable :: x(:)
  type(t_perm) :: perm
  type(t_perm), allocatable :: ptab(:)
  type(t_dperm) :: dperm, dper2

   ! build the permutation table
   do d=1,4
     call perm_table(ptab,d)
     write(*,*) "Permutation table for d=",d,":"
     do i=1,size(ptab)
       write(*,*) "  permutation ",i
       write(*,*) "    %d = ",ptab(i)%d
       write(*,*) "    %p = ",ptab(i)%p
       write(*,*) "    %i = ",ptab(i)%i
     enddo
     write(*,*) "  idx(permutation table) = ",idx(ptab), &
         " (should be idx.eq.[1,...,d] )"
     deallocate(ptab)
   enddo

   ! test functions working on t_dperm objects
   write(*,*) "Data permutation tests"
   do d=1,10
     allocate(v(d),x(d)); call random_number(x); v = int(100.0*x)
     ! generate some repetitions
     if(d.ge.5) v(d/2) = v(d)
     if(d.ge.7) v(d/3) = v(d)
     dperm = dperm_reduce(v)
     write(*,*) "  data  = ",v
     write(*,*) "  dperm%pi%d = ",dperm%pi%d
     write(*,*) "  dperm%pi%p = ",dperm%pi%p
     write(*,*) "  dperm%pi%i = ",dperm%pi%i
     write(*,*) "  dperm%x    = ",dperm%x
     write(*,*) "  perm error = ",maxval(abs(dperm%x-v(dperm%pi%i)))
     deallocate(v,x)
   enddo

   ! test the operators
   write(*,*) "Operator tests"
   deallocate(dperm%x)
   do d=1,10
     allocate(v(d),x(d))
     call random_number(x); v = int(10.0*x)-5
     dperm%pi%d = d; allocate(dperm%x(d)); dperm%x = v
     if(d.ne.3) & ! one case of data1.eq.data2
       call random_number(x); v = int(4.0*x)-2
     dper2%pi%d = d; allocate(dper2%x(d)); dper2%x = v
     write(*,*) "  data1 = ",dperm%x
     write(*,*) "  data2 = ",dper2%x
     write(*,*) "    data1.eq.data2 = ",dperm.eq.dper2
     write(*,*) "    data1.ne.data2 = ",dperm.ne.dper2
     write(*,*) "    data1.lt.data2 = ",dperm.lt.dper2
     write(*,*) "    data1.le.data2 = ",dperm.le.dper2
     write(*,*) "    data1.gt.data2 = ",dperm.gt.dper2
     write(*,*) "    data1.ge.data2 = ",dperm.ge.dper2
     deallocate(v,x,dperm%x,dper2%x)
   enddo
   do d=1,7
     ! permutation operators
     call perm_table(ptab,d)
     allocate(v(2),x(2))
     call random_number(x); v = int(real(fact(d))*x)
     where(v.eq.0) v=1
     write(*,*) "  perm1 = ",ptab(v(1))%i
     write(*,*) "  perm2 = ",ptab(v(2))%i
     perm = ptab(v(1))*ptab(v(2))
     write(*,*) "    perm1*perm2 = ",perm%p,perm%i
     perm = ptab(v(1))**2
     write(*,*) "    perm1**2    = ",perm%p,perm%i
     perm = ptab(v(1))**(-1)
     write(*,*) "    perm1**(-1) = ",perm%p,perm%i
     perm = ptab(v(1))*(ptab(v(1))**(-1))
     write(*,*) "    perm1*perm1**(-1) = ",perm%p,perm%i
     deallocate(v,x)
   enddo
   deallocate(ptab)

   ! test the combinations
   call comb_table(combs,(/ 5 , 3 , 8 , 12 , 1 , 0 /),3)
   write(*,*) "Combination table:"
   write(*,*) "  shape: ",shape(combs)
   write(*,'(3i4)') combs

 end subroutine test_module

!-----------------------------------------------------------------------

end module mod_perms

